Drift and trapping in biased diffusion on disordered lattices 
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o\ : 

i We reexamine the theory of transition from drift to no-drift in biased diffusion on percolation 

■ networks. We argue that for the bias field B equal to the critical value B c , the average velocity 

- .. ' at large times t decreases to zero as l/log(t). For B < B c , the time required to reach the steady- 

state velocity diverges as exp(const/|B c — B\). We propose an extrapolation form that describes 
the behavior of average velocity as a function of time at intermediate time scales. This form is 
found to have a very good agreement with the results of extensive Monte Carlo simulations on a 
^"■^ , 3-dimensional site-percolation network and moderate bias. 

i—i, Diffusion in disordered lattices and in particular at their percolation threshold is an old subject [0. If there is 
' a preferred direction of motion due to some external imposed field, we have biased diffusion. A small value of the 
, external bias gives rise to a mean displacement in the direction of field that increases linearly with time, and the mean 
velocity tends to a constant. This asymptotic value of the mean velocity is proportional to the field for small fields. 
I , It was pointed out by Bottger and Bryksin Q that in disordered media, because of trapping in dead-end branches, 
the mean velocity would be a nonmonotonic function of the bias. In it was argued that for random walkers with 
no mutual interactions, the asymptotic mean velocity actually becomes zero for a finite value of bias B c , and that the 
mean displacement increases as a sublinear power of time for B > B c . These theoretical arguments were supported by 
exact calulation of the mean drift velocity on a random comb |5| . A similar argument was used to show the existence 
of drift to no drift transition in the case of non- interacting particles when the bias field direction is not constant, but 
depends on the random geometry || . It was argued in [Q that this sharp transition from drift to no drift disappears 
if a repulsive interaction between diffusing particles is taken into account. 

However, clear numerical verification of these theoretical predictions has not been possible so far. Earlier numerical 
simulations failed to see clear evidence of a sharp transition from drift to no drift as a function of the bias |^j . For large 
bias, log-periodic oscillations in time mask the possible transition There are large sample to sample fluctuations, 
and it appears that times much greater than feasible are needed to see the asymptotic time regime for B near B c . 
, This has motivated our reexamination of this question in this paper. 

In this paper, we argue that the time needed to see the asymptotic behavior predicted in J3|,|j diverges as 
exp(const/|B c — B\), as B approaches B c from below. This kind of sharp increase of relaxation times, often en- 
countered in systems with quenched disorder (the familiar Vogel-Fulcher law in glasses), implies that in order to 
compare the theoretical predictions to simulations results, one must take into account corrections to the true asymp- 
qq 1 totic behavior. We argue that at B — B c , the mean velocity decreases as 1/logi for large times t. We propose a 
^\ , simple extrapolation form that incorporates this behavior, and agrees with the expected asymptotic behavior of mean 
1 velocity for large times. We test this form by comparing to results of Monte Carlo simulations on a 3-d simple cubic 
lattice. We find excellent agreement between the two. 

We start with a precise definition of the problem, and of the method used in computer simulations. We consider a 
site percolation problem on an L x L x L simple cubic lattice with periodic boundary conditions with concentration 
p — 1/2 of sites colored white, and the rest red. This value of p is much larger than the critical percolation threshold 
p c = 0.31160, and for this value of p, more that 97.8% of the white sites belong to the infinite white cluster. At time 
O ■ t = 0, a group of N random walkers are placed on randomly selected white sites of the lattice. At each (discrete) 
time step, each walker attempts to move to a nearest neighbor site, and actually moves there if the site is white. The 
nearest neighbor is chosen to be the neighbour along the positive x-direction with bias probability B, and to any of 
the six neighbours with probability 1 — 5; thus B+ (1 — B)/6 = (1 + 5B)/6 is the probability to move into the positive 
j3 ■ x direction, and (1 — B)/6 is the one for the negative x direction. As in the earlier simulations ||, the status of lattice 
sites was stored in single bits; thus L = 964 could be chosen for most of our simulations. The total displacement of a 
walker is calculated by adding up all her single-step displacements. This is then averaged over the N walkers. In our 
simulations, N varies from 80,000 for the shorter time runs, to only 384 for the longest simulations. 

For weak bias B, the situation is quite straightforward. The average velocity of the walkers quickly settles down 
to a constant value which for small bias increases linearly with bias. In fig 1, we have plotted the limiting velocity as 
a function of B, for B lying between 0.15 and 0.47. For larger B, the velocity decreases very slowly with time. We 
can not be sure of the true asymptotic velocity. The plotted mean velocities are determined by looking at the rms 
displacements for 0.9 x 10 9 <i<2.1x 10 9 . Also plotted in the figure are the limiting values obtained by using the 
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extrapolation formula (see below) for B > .4. For lower values of B, these match very well. But for B = .45 and 
above, clearly much longer times would be needed to reach the asymptotic values. 

In the presence of strong bias, the diffusive motion of walkers in the disordered lattice is slowed down because often 
the walker gets into a cul-de-sac region, and it takes a long time to get out of them. For a dead-end of depth £ the 
trapping time varies as X e , where A is the ratio of transition probabilities along and against the field. In our case 
A = (1 + 5B)/(l — B). In addition, according to Refs. the density of traps of depth £ decreases exponentially 
with I, say Prob(£) varies as exp(— At), where A is a p— dependent constant. The average trapping time < t > at a 
backbone site is given by MM 



^Prob(^)A £ (1) 



For B > B c , this summation diverges, and the asymptotic drift velocity , which varies inversely as the trapping time, 
is zero. It is easy to see that the critical value B c is determined by the equation 

(1 + 5£ C )/(1-B c )exp(-A) = l (2) 

For B < B c , the average trapping time is finite, and the asymptotic drift velocity is finite, and decreases to zero 
as B tends to B c from below. For B > B c , the walker moves a distance R of order l/Prob(£) before it meets its first 
encountered trap of depth £. The time t to escape from it is of the order of A^. Thus from eq(2) we have Q 

R oc exp(Af) = exp(yllnVlnA) = exp(lntln A c /ln A) = t lnX ^ /lnX 

Thus the average distance moved increases as t 1 ~ x where 

1 log[(l + 5£? c )/(l - B c )] 



log[(l + 5B)/(l-B)] 



(3) 



If p = p c , the critical threshold for percolation, then the constant A, and hence using eq. (2), B c is zero. In this 
case, x = 1, and the mean displacement grows only as a power of logt. This has been observed in simulations [ p"o| . 

Consider now the motion of a large number of walkers at B = B c at large times t. Let us estimate the average 
trapping time felt by one walker in this time. This is approximately given by Eq. (1), except that the summation 
is cut off at a value £ m ax\ which is determined by the condition that A^ max is order t. This is because in this time 
the walker is unlikely to have encountered a deeper trap. As each term of this summation is roughly equal this 
implies that average trapping time up to time t increases as £ m ax ~ log(i)- Using the fact that average velocity is just 
the inverse of the average trapping time, we see that 

v(t) ~ 1/logi, for B = B c . (4) 

For B = B c — e, the summation in Eq. (1) is finite, but diverges as 1/e for small e. This implies that v(t) tends to 
a constant proportional to (B c — B) for B < B c for large t. It varies as 1/ logi for B = B c and large t, and varies as 
t~ x for B > B c , with x determined by Eq. (3). A simple extrapolation form which incorporates all these behaviors is 

v = Kx/[(t/t o r - 1] (5) 

where if is a constant, independent of t, but weakly dependent on the bias field B. If x is positive, v decreases as 
t~ x for large t. For x — it varies as 1/ logt, and for x < 0, it tends to a finite limit as t — > oo. For B = B c — e, with 
e small, initially, up to some relaxation time T, the velocity will decrease roughly as 1/logi, but for t > T it reaches 
the constant asymptotic value proportional to e. Matching these two values of velocity at t = T, we get 

T ~ exp(const/e) (6) 

Thus we see that the relaxation time of the system for B near B c diverges according to the Vogel-Fulcher form, as 
was claimed in the introduction. 

Eq.(5) suggests a simple way to analyse the simulation data to find B c . We plot l/v versus \ogt for various values of 
the bias field B. The plot is linear right at the possible transition point, shows a decreasing slope, eventually settling 
to a finite constant value for smaller B, and shows slope increasing with time for larger B. This is shown in Fig. 2. 
We find that for B = 0.53, the graph is fairly linear. The equation of the best fit straight line is 
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1/v = 11.2 log(t/30) 



(7) 



Comparing with the extrapolation formula, we see that for x —> 0, the right hand side reduces to Kj log(t/to)- This 
fixes the parameters B c ~ 0.53, and to — 30. We have shown by continuous curves theoretical fits to the data in Fig 
2, for other values of B, using Eq. (3) to determine x for a given choice of B. The only unknown parameter K is 
expected to depend on B weakly. We determine its value by selecting the best fit to the data. The best fit value of 
K is found to be 0.0678 for £=0.40, which increases to 0.0893 for B=.53 and 0.0927 for B=0.60. We see that by a 
suitable choice of K, the extrapolation formula (5) provides a very good fit to the data in the entire range of data 
t > 10 3 . For t < 10 3 , there are significant corrections due to short term transients, not taken into account in the 
extrapolation form (5). 

To better test Eq. (4), we have run the simulations precisely for B = B c for much longer times. The results are 
shown in Fig 3. We see that no significant deviations from the linear behavior are seen for t up to > 10 10 . Note the 
large fluctuations for large time data, which is averaged over only a small number of walkers. For shorter times these 
results were confirmed by P. Grassberger (private communication) using a very different algorithm. 

We see in Fig. 2 that with two global parameters ( B c and to), and only one free parameter K for each bias value, 
we can fit nicely the data for t > 10 3 for all the bias fields 0.4 < B < 0.6 to the functional form Eq. (5), making 
the above conclusions based on it more trustworthy. Thus our data are in good agreement with the theoretically 
predicted functional form, and provide the first direct observation in simulations of a sharp transition between drift 
and no drift. We expect a similar behavior to be present in the case of 'topological bias' studied in ||. 
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FIG. 1. The mean velocity at large times plotted against the bias field B. The observed values in simulations at times of 
the order 10 9 are shown by diamonds. The extrapolated limiting values are shown by crosses. Times are measured in units of 
diffusion attempts per walker, and distances in units of the lattice constant. 




FIG. 2. Reciprocal velocity versus log(time) showing the transition from concave to convex curvature at B c ~ 0.53 for 
intermediate times. 




FIG. 3. Reciprocal velocity versus log(time) for B = B c . Different symbols show data averaged over different number N of 
walkers. N = 80,000 (o), 64,000 (+), 20000 (squares), 1024 (x), 384 (A). 
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